Results for Supplemental Material for Hitchcock & Frank 2023
Peter Frank Hitchcock - LNCC Lab, Brown University
2023-06-23
sapply(c(
"rjson",
"data.table",
"dplyr",
"ggplot2",
"stringr",
"purrr",
"foreach",
"doParallel",
"patchwork",
"lme4",
"lmerTest",
"testit",
"latex2exp",
"ggpubr"
),
require, character=TRUE)## Loading required package: rjson## Loading required package: data.table## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: ggplot2## Loading required package: stringr## Loading required package: purrr##
## Attaching package: 'purrr'## The following object is masked from 'package:data.table':
##
## transpose## Loading required package: foreach##
## Attaching package: 'foreach'## The following objects are masked from 'package:purrr':
##
## accumulate, when## Loading required package: doParallel## Loading required package: iterators## Loading required package: parallel## Loading required package: patchwork## Loading required package: lme4## Loading required package: Matrix## Loading required package: lmerTest##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step## Loading required package: testit## Loading required package: latex2exp## Loading required package: ggpubr## rjson data.table dplyr ggplot2 stringr purrr foreach
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## doParallel patchwork lme4 lmerTest testit latex2exp ggpubr
## TRUE TRUE TRUE TRUE TRUE TRUE TRUEsf <- function() sapply(paste0("./Functions/", list.files("./Functions/", recursive=TRUE)), source) # Source all fxs
sf()## ./Functions/General/FindDelay.R ./Functions/General/FindDelayAlt.R
## value ? ?
## visible FALSE FALSE
## ./Functions/General/Utils.R ./Functions/Model/Models/RunMBayesLearner.R
## value ? ?
## visible FALSE FALSE
## ./Functions/Model/Models/RunMBayesLearnerDiffBeta.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayes.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayesAlt.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayes.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayesAlt.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearner.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayTo0Inits.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPrior.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffBeta.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffLR.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPrior.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorESAndEpsFixed.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorNoCK.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/AltPlotSimTrainingCurves.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/aModelUtils.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotAllWorstTestOptionsSimVsEmp.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotSimEmpTest.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotSimEmpTestNoSim.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PrepDfForModel.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PrepDfForModelPROpt.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/RecodeDfs.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/SimplePlotSimTrainingDelay.R
## value ?
## visible FALSE
## ./Functions/Model/OptFxs/RunOptTrainSIT.R
## value ?
## visible FALSE
## ./Functions/Model/OptFxs/RunOptTrainSITPR.R
## value ?
## visible FALSE
## ./Functions/Model/SimFxs/RunSimTrainSIT.R
## value ?
## visible FALSE
## ./Functions/Model/SimFxs/RunSimTrainSITForPR.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/aModelFunctions.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesMean.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesStd.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesVar.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/DecayChoiceKernel.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/DecayQVals.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/UpdateABMatrix.R
## value ?
## visible FALSEDefPlotPars()
registerDoParallel(cores=round(detectCores()*2/3))Load data from studies 1 and 2
s1_train <- read.csv("../data/cleaned_files/s1_train_with_delay_deident.csv")
s1_sit <- read.csv("../data/cleaned_files/s1_SIT_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_train <- read.csv("../data/cleaned_files/s2_train_with_delay_deident.csv")
s2_sit <- read.csv("../data/cleaned_files/s2_sit_deident_corrected_names.csv")
s1_pst <- read.csv("../data/cleaned_files/s1_PST_deident.csv") %>% rename(ID=deident_ID)
s2_pst <- read.csv("../data/cleaned_files/s2_PST_deident.csv") %>% rename(ID=deident_ID)Harmonize variables and create some separate vars
# Study 2 harmonize
s2_sit[s2_sit$condition=="cogn", "condition"] <- "cognitive"
s2_train[s2_train$trial_within_condition <= 20, "block"] <- 1
s2_train[s2_train$trial_within_condition > 20, "block"] <- 2
s1_sit$probability <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 1)))
s1_sit$valence <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 2)))
s2_sit$probability <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 1)))
s2_sit$valence <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 2)))Define paths and functions
# MLE results
allp <- "../model_res/opts_mle_paper_final/all/"
# Empirical Bayes results
bp <- "../model_res/opts_mle_paper_final/best/"
# Simulations
sp <- "../model_res/sims_clean/sims_from_mle/"
# Read model function
rm <- function(path, model_str) read.csv(paste0(path, model_str))Plotting delay
s1_train_delay_summs <- s1_train %>% filter(!is.na(delay_trunc)) %>% group_by(condition, delay_trunc, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition', 'delay_trunc'. You can
## override using the `.groups` argument.s1_tr_delay_summs <-
Rmisc::summarySEwithin(s1_train_delay_summs,
measurevar = "m",
withinvars = c("condition", "delay_trunc"),
idvar = "ID")## Automatically converting the following non-factors to factors: condition, delay_truncs2_train_delay_summs <- s2_train %>% filter(!is.na(delay_trunc)) %>% group_by(condition, delay_trunc, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition', 'delay_trunc'. You can
## override using the `.groups` argument.s2_tr_delay_summs <-
Rmisc::summarySEwithin(s2_train_delay_summs,
measurevar = "m",
withinvars = c("condition", "delay_trunc"),
idvar = "ID")## Automatically converting the following non-factors to factors: condition, delay_trunca <- ggplot(s1_tr_delay_summs, aes(x=delay_trunc, y=m, group=condition, color=condition)) +
geom_line(size=2, position = position_dodge(width = .2)) +
geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3, color="gray57") +
geom_hline(yintercept = .5, size=2, color="black") +
geom_errorbar(aes(x=delay_trunc, y=m, ymin=m-se, ymax=m+se, group=condition),
position = position_dodge(width = .2), color="black", size=1.5, width=.15) +
geom_point(size=6, position = position_dodge(width = .2), pch=21, fill="gray57", alpha=.95) +
ga + ap + ft + tol + ylim(.48, .82) + tp +
scale_color_manual(values=c("purple", "orange"),
labels=c("cognitive", "overt")) + ylab("Proportion correct") + xlab("Delay") +
ggtitle("Experiment 1")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.b <- ggplot(s2_tr_delay_summs, aes(x=delay_trunc, y=m, group=condition, color=condition)) +
geom_line(size=2, position = position_dodge(width = .2)) +
geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3, color="gray57") +
geom_hline(yintercept = .5, size=2, color="black") +
geom_errorbar(aes(x=delay_trunc, y=m, ymin=m-se, ymax=m+se, group=condition),
position = position_dodge(width = .2), color="black", size=1.5, width=.15) +
geom_point(size=6, position = position_dodge(width = .2), pch=21, fill="gray57", alpha=.95) +
ga + ap + ft + lp + ylim(.48, .82) + tp +
scale_color_manual(values=c("purple", "orange"),
labels=c("cognitive", "overt")) + ylab("") + xlab("Delay") +
ggtitle("Experiment 2")delay_comb_plot <- a+bdelay_comb_plot
#ggsave("../paper/figs/fig_parts/delay_plot.png", delay_comb_plot, height=7, width=10, dpi=700)Supplemental material model validation (although learning curves exp 2 plot is in key_results)
6/5/23 — reran sims to mkae sure completely updated/bug fixed version of par results used
s1_train_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDiffDecayToPessPrior57088.csv")
s2_train_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_2_train_SIT__train_RunMQLearnerDiffDecayToPessPrior28414.csv")
s1_test_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior57088.csv")
s2_test_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior28414.csv")Capture of delay effects
out_a <- SimplePlotSimTrainingDelay(emp_df=s1_train, sim_df=s1_train_sim_m1_eb)## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'iter'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.## Warning: Duplicated aesthetics after name standardisation: colourout_b <- SimplePlotSimTrainingDelay(emp_df=s2_train, sim_df=s2_train_sim_m1_eb) ## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'iter'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.## Warning: Duplicated aesthetics after name standardisation: colourout_a## Warning: Duplicated aesthetics after name standardisation: colour
out_b## Warning: Duplicated aesthetics after name standardisation: colour
# ggsave("../paper/figs/fig_parts/simvsemp_delay_plot_1.png", out_a, height=8, width=12, dpi=700)
# ggsave("../paper/figs/fig_parts/simvsemp_delay_plot_2.png", out_b, height=8, width=12, dpi=700)Worst-option Proportions in Test Phase and comparison to same model without CK
worst_plots_a_ck <- PlotAllWorstTestOptionsSimVsEmp(s1_sit, s1_test_sim_m1_eb, m_string="Simulations \nwith Choice Kernel")w_a_ck <- worst_plots_a_ck[[1]]+worst_plots_a_ck[[2]]w_a_ck 
worst_plots_b_ck <- PlotAllWorstTestOptionsSimVsEmp(s2_sit, s2_test_sim_m1_eb, m_string="Simulations \nwith Choice Kernel")w_b_ck <- worst_plots_b_ck[[1]]+worst_plots_b_ck[[2]]w_b_ck## Warning: Removed 1 rows containing missing values (`geom_point()`).
# ggsave("../paper/figs/fig_parts/CK_simvsemp_worst_plot_1.png", w_a_ck, height=6, width=15, dpi=700)
# ggsave("../paper/figs/fig_parts/CK_simvsemp_worst_plot_2.png", w_b_ck, height=6, width=15, dpi=700)Comparative worst plots
s1_test_sim_m1_eb_no_ck <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPriorNoCK16278.csv")
s2_test_sim_m1_eb_no_ck <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPriorNoCK31508.csv")worst_plots_nck_a <- PlotAllWorstTestOptionsSimVsEmp(s1_sit, s1_test_sim_m1_eb_no_ck, m_string="Simulations \nNo Choice Kernel")w_a_nck <- worst_plots_nck_a[[1]]+worst_plots_nck_a[[2]]w_a_nck
worst_plots_nck_b <- PlotAllWorstTestOptionsSimVsEmp(s2_sit, s2_test_sim_m1_eb_no_ck, m_string="Simulations \nNo Choice Kernel")Captures the mis-spec
w_b_nck <- worst_plots_nck_b[[1]]+worst_plots_nck_b[[2]]w_b_nck
# ggsave("../paper/figs/fig_parts/NOCK_simvsemp_worst_plot_1.png", w_a_nck, height=6, width=15, dpi=700)
# ggsave("../paper/figs/fig_parts/NOCKsimvsemp_worst_plot_2.png", w_b_nck, height=6, width=15, dpi=700)Parameter recovery
m1_par_recov <-
read.csv("../model_res/par_recov_clean/pr_opts_mle/best_pr/par_recov_BEST_EMPIRICAL_BAYES_mv_gaussstudy_1_train_SIT_RunMQLearnerDiffDecayToPessPrior_76217.csv") # Confirmed created 3/28 after the 3/23 bug fix on decay eb_recov_pars_m1 <- m1_par_recov[grep("EB_recovered", names(m1_par_recov))]
simmed_pars_m1 <- m1_par_recov[grep("simmed", names(m1_par_recov))]Phi difference
this_eb_pr_df <-
data.frame("simulated"=simmed_pars_m1$cog_phi_simmed-simmed_pars_m1$overt_phi_simmed,
"EB_recovered"=eb_recov_pars_m1$cog_phi_EB_recovered-eb_recov_pars_m1$overt_phi_EB_recovered)
phi_diff <- ggplot(this_eb_pr_df, aes(x=simulated, y=EB_recovered)) +
geom_line(aes(x=simulated, y=simulated), size=4) +
geom_point(size=8, pch=21, fill="gray57", alpha=.8) +
stat_cor(method="spearman", size=6, label.y=.9) +
ggtitle(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
ga + ap + tp +
ylab("Recovered") + xlab("Simulated")
phi_diff
#ggsave("../paper/figs/fig_parts/pr/pr_phi_diff.png", phi_diff, height=4, width=8, dpi=700)Percent correctly recovering the sign difference
this_eb_pr_df$sim_sign <- if_else(this_eb_pr_df$simulated >= 0, 1, 0)
this_eb_pr_df$eb_sign <- if_else(this_eb_pr_df$EB_recovered >= 0, 1, 0)
# this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign
# (this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign)*1
ta <- table((this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign)*1)
ta[2]/sum(ta)## 1
## 0.7142857plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$epsilon_simmed, eb_recov_pars_m1$epsilon_EB_recovered,
"$\\epsilon", stat_pos = .32)
plot
# ggsave("../paper/figs/fig_parts/pr/eps.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$cog_phi_simmed, eb_recov_pars_m1$cog_phi_EB_recovered,
"$\\phi^{Cog}")
plot
# ggsave("../paper/figs/fig_parts/pr/phi_cog.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$overt_phi_simmed, eb_recov_pars_m1$overt_phi_EB_recovered,
"$\\phi^{Overt}")
plot
# ggsave("../paper/figs/fig_parts/pr/phi_overt.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$choice_LR_simmed, eb_recov_pars_m1$choice_LR_EB_recovered,
"$\\alpha_{CK}")
plot
# ggsave("../paper/figs/fig_parts/pr/ck.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$explor_scalar_simmed, eb_recov_pars_m1$explor_scalar_EB_recovered,
"ES")
plot
# ggsave("../paper/figs/fig_parts/pr/es.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$q_LR_simmed, eb_recov_pars_m1$q_LR_EB_recovered,
"$\\alpha_{Q}")
plot
# ggsave("../paper/figs/fig_parts/pr/qlr.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$beta_simmed, eb_recov_pars_m1$beta_EB_recovered,
"$\\beta", stat_pos = 25)
plot
Simulation with a high difference in phi
Load model fits
m1_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior17352.csv")
m1_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior58319.csv")
m1_study1_eb <- rbind(m1_study1_eb_v1, m1_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m1_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior12123.csv")
m1_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior74336.csv")
m1_study2_eb <- rbind(m1_study2_eb_v1, m1_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))Find a participant with relatively high diff in decay
m1_study1_eb$cog_phi[25]## [1] 0.6817201m1_study1_eb$overt_phi[25]## [1] 0.1300788m1_study1_eb[25, ]## # A tibble: 1 × 12
## # Groups: ID [1]
## X epsilon q_LR cog_phi overt_phi choice_LR explor_scalar beta
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 25 0.0365 0.732 0.682 0.130 0.0425 0.220 4.97
## # … with 4 more variables: convergence <int>, nll <dbl>, AIC <dbl>, ID <int>This did correspond to a big difference in performance empirically
s1_train %>% filter(ID==25) %>% group_by(valence_and_probability, condition) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'valence_and_probability'. You can override
## using the `.groups` argument.## # A tibble: 8 × 3
## # Groups: valence_and_probability [4]
## valence_and_probability condition m
## <chr> <chr> <dbl>
## 1 40-10_punishment cognitive 0.6
## 2 40-10_punishment overt 0.95
## 3 40-10_reward cognitive 0.9
## 4 40-10_reward overt 1
## 5 90-10_punishment cognitive 0.525
## 6 90-10_punishment overt 0.775
## 7 90-10_reward cognitive 0.05
## 8 90-10_reward overt 1s1_sit %>% filter(ID==25) %>% group_by(valence_and_probability, condition) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'valence_and_probability'. You can override
## using the `.groups` argument.## # A tibble: 8 × 3
## # Groups: valence_and_probability [4]
## valence_and_probability condition m
## <chr> <chr> <dbl>
## 1 40-10_punishment cognitive 0.75
## 2 40-10_punishment overt 1
## 3 40-10_reward cognitive 1
## 4 40-10_reward overt 1
## 5 90-10_punishment cognitive 1
## 6 90-10_punishment overt 1
## 7 90-10_reward cognitive 0
## 8 90-10_reward overt 1Summarize the correct and incorrect q-values
corr_q_vals_summ <-
s1_train_sim_m1_eb %>% filter(ID==25) %>%
group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_correct_Q_vals))## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.corr_q_vals_summ$valence <- factor(corr_q_vals_summ$valence, levels=c("reward", "punishment"))
corr_q_vals_summ[corr_q_vals_summ$cond=="cog", "cond"] <- "cognitive"incorr_q_vals_summ <-
s1_train_sim_m1_eb %>% filter(ID==25) %>%
group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_incorrect_Q_vals))## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.incorr_q_vals_summ$valence <- factor(incorr_q_vals_summ$valence, levels=c("reward", "punishment"))
incorr_q_vals_summ[incorr_q_vals_summ$cond=="cog", "cond"] <- "cognitive"Summarize the simulated performance during learning
sim_perf_summ <-
s1_train_sim_m1_eb %>% filter(ID==25) %>%
group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_corrs))## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.sim_perf_summ$valence <- factor(sim_perf_summ$valence, levels=c("reward", "punishment"))
sim_perf_summ[sim_perf_summ$cond=="cog", "cond"] <- "cognitive"Take the difference in q-values
assert(all(corr_q_vals_summ$trial_within_condition==incorr_q_vals_summ$trial_within_condition))
assert(all(corr_q_vals_summ$valence==incorr_q_vals_summ$valence))
assert(all(corr_q_vals_summ$probability==incorr_q_vals_summ$probability))
corr_q_vals_summ$diff <- corr_q_vals_summ$m-incorr_q_vals_summ$ma <-
ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),
aes(x=trial_within_condition, y=m, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ cond) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.9, color="blue", size=1.5, linetype="longdash") +
geom_hline(yintercept=-.1, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("Q-value") +
ggtitle("Correct") +
xlab("Stim iteration") + ylim(-1, 1)b <-
ggplot(incorr_q_vals_summ %>% filter(probability=="90-10"),
aes(x=trial_within_condition, y=m, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ cond) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.1, color="blue", size=1.5, linetype="longdash") +
geom_hline(yintercept=-.9, color="red", size=1.5, linetype="longdash") +
ga + ap + tp + ft + tol +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ggtitle("Incorrect") +
ylab("") +
xlab("Stim iteration") + ylim(-1, 1)c <-
ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),
aes(x=trial_within_condition, y=diff, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ cond) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="blue", size=1.5, linetype="longdash") +
geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Difference") +
xlab("Stim iteration") + ylim(-1, 1)
# For legend
# ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),
# aes(x=trial_within_condition, y=diff, color=valence)) +
# geom_hline(yintercept=0, color="black", size=2) +
# geom_line(size=3) + facet_wrap( ~ cond) +
# scale_color_manual(values=c("blue", "red")) +
# geom_hline(yintercept=.8, color="blue", size=1.5, linetype="longdash") +
# geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
# ga + ap + ft +
# theme(legend.text = element_text(size = 30),
# legend.title = element_blank(),
# legend.key.size = unit(2.5, 'lines')) +
# theme(plot.title = element_text(size = 30, hjust = .5)) +
# ylab("") +
# ggtitle("Difference") +
# xlab("Stim iteration") + ylim(-1, 1)d <- ggplot(sim_perf_summ %>% filter(probability=="90-10"),
aes(x=trial_within_condition, y=m, color=valence)) +
geom_hline(yintercept=0.5, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ cond) +
scale_color_manual(values=c("blue", "red")) +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("Proportion correct") +
ggtitle("Learning") +
xlab("Stim iteration") sim_test_perf_summ <-
s1_test_sim_m1_eb %>% filter(ID==25) %>%
group_by(cond, probability, valence) %>% summarize(m=mean(sit_sim_corrs))## `summarise()` has grouped output by 'cond', 'probability'. You can override
## using the `.groups` argument.sim_test_perf_summ$valence <- factor(sim_test_perf_summ$valence, levels=c("reward", "punishment"))
sim_test_perf_summ[sim_test_perf_summ$cond=="cog", "cond"] <- "cognitive"#sim_test_perf_summ %>% filter(probability=="90-10")
e <- ggplot(sim_test_perf_summ %>% filter(probability=="90-10"),
aes(x=valence, y=m, color=valence)) +
geom_hline(yintercept=0.5, color="black", size=2) +
geom_bar(stat="identity", fill="white", size=3) +
facet_wrap(~cond) +
scale_color_manual(values=c("blue", "red")) +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("Proportion correct") +
ggtitle("Test") +
xlab("") + theme(axis.ticks.x=element_blank(), axis.text.x=element_blank())#+ ylim(.48, 1) +q_vals_comb <- a + b + cq_vals_comb
sim_perf_comb <- d + esim_perf_comb
#ggsave("../paper/figs/fig_parts/q_vals_plot.png", q_vals_comb, height=5, width=16, dpi=700)#ggsave("../paper/figs/fig_parts/sim_perf_comb.png", sim_perf_comb, height=6, width=14, dpi=700)Supplementary Figure 7
Plotting Q-values by valence for m3 vs. m4 (m6 vs. m7) to understand how the introduction of pessimistic priors affects performance
# Pessimistic
m3_sims_s1 <-
rm(sp, "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDecayToPessPrior82512.csv")
# Naive
m4_sims_s1 <-
rm(sp, "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDecayTo0Inits36394.csv")M6: Q nothing varied in paper
pess_priors_qv_summs <- m3_sims_s1 %>% filter(probability=="90-10") %>%
group_by(valence, trial_within_condition) %>%
summarize(m_corr=mean(sim_correct_Q_vals), m_incorr=mean(sim_incorrect_Q_vals)) %>%
mutate(diff=m_corr-m_incorr)## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.pess_priors_qv_summs$valence <- factor(pess_priors_qv_summs$valence, levels=c("reward", "punishment"))M7: Q nothing varied, decay to 0
naive_priors_qv_summs <- m4_sims_s1 %>% filter(probability=="90-10") %>%
group_by(valence, trial_within_condition) %>%
summarize(m_corr=mean(sim_correct_Q_vals), m_incorr=mean(sim_incorrect_Q_vals)) %>%
mutate(diff=m_corr-m_incorr)## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.naive_priors_qv_summs$valence <- factor(naive_priors_qv_summs$valence, levels=c("reward", "punishment"))Break down into correct/incorrect plots
ggplot(pess_priors_qv_summs,
aes(x=trial_within_condition, y=m_corr, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
# geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Correct") +
xlab("Stim iteration") + ylim(-1, 1)
ggplot(naive_priors_qv_summs,
aes(x=trial_within_condition, y=m_corr, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
# geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Correct") +
xlab("Stim iteration") + ylim(-1, 1)
ggplot(pess_priors_qv_summs,
aes(x=trial_within_condition, y=m_incorr, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
#geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
# geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Incorrect") +
xlab("Stim iteration") + ylim(-1, 1)
ggplot(naive_priors_qv_summs,
aes(x=trial_within_condition, y=m_incorr, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
# geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Incorrect") +
xlab("Stim iteration") + ylim(-1, 1)
pess_plot <- ggplot(pess_priors_qv_summs,
aes(x=trial_within_condition, y=diff, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("Q-value Difference: \n Correct vs. Incorrect") +
ggtitle("Pessimistic Priors") +
xlab("Stim iteration") + ylim(0, 1)naive_plot <- ggplot(naive_priors_qv_summs,
aes(x=trial_within_condition, y=diff, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Naive Priors (All Zero)") +
xlab("Stim iteration") + ylim(0, 1)simple_qv_plots_comb <- pess_plot + naive_plotsimple_qv_plots_comb 
No powerpoint — loaded directly into supplemental
#ggsave("../paper/figs/fig_parts/simple_qv_plots_comb.png", simple_qv_plots_comb, height=7, width=14, dpi=700)PST/Stimulus valuation phase analyses
Used rew history because Q-values are for specific Q(s,a)s whereas selection in the PST phase is between stimuli ie. Q(s, :)
# Take just the trials within condition
within_cond_pst_s1 <- s1_pst %>% filter(test_condition == "cogn_cogn" | test_condition == "overt_overt")
within_cond_pst_s1$test_condition <- factor(within_cond_pst_s1$test_condition, levels=c("overt_overt", "cogn_cogn"))
contrasts(within_cond_pst_s1$test_condition) <- c(-.5, .5)
head(within_cond_pst_s1$test_condition)## [1] overt_overt overt_overt overt_overt overt_overt overt_overt overt_overt
## attr(,"contrasts")
## [,1]
## overt_overt -0.5
## cogn_cogn 0.5
## Levels: overt_overt cogn_cogn# Both singular
# summary(m0_s1 <- glmer(resp_num ~ scale(left_min_right) + ( scale(left_min_right)|ID),
# data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))
#
# summary(m0_s1 <- glmer(resp_num ~ scale(left_min_right)*test_condition + (scale(left_min_right) |ID),
# data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))Using model without random intercepts because not singular
summary(m0_s1 <- glmer(resp_num ~ scale(left_min_right) + (0 + scale(left_min_right)|ID),
data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: resp_num ~ scale(left_min_right) + (0 + scale(left_min_right) |
## ID)
## Data: within_cond_pst_s1
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 6456.5 6476.6 -3225.2 6450.5 5983
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3395 -0.7468 0.1149 0.7468 4.5130
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID scale(left_min_right) 1.538 1.24
## Number of obs: 5986, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.11210 0.03073 3.648 0.000264 ***
## scale(left_min_right) 1.41202 0.12084 11.685 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## scl(lft_m_) 0.015summary(m1_s1 <- glmer(resp_num ~ scale(left_min_right)*test_condition + (0 + scale(left_min_right) |ID),
data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## resp_num ~ scale(left_min_right) * test_condition + (0 + scale(left_min_right) |
## ID)
## Data: within_cond_pst_s1
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 6431.2 6464.7 -3210.6 6421.2 5981
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.1380 -0.7457 0.1136 0.7448 4.1139
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID scale(left_min_right) 1.555 1.247
## Number of obs: 5986, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.11079 0.03085 3.592 0.000328 ***
## scale(left_min_right) 1.42529 0.12150 11.731 < 2e-16 ***
## test_condition1 0.09218 0.06167 1.495 0.134955
## scale(left_min_right):test_condition1 -0.36515 0.07029 -5.195 2.05e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sc(__) tst_c1
## scl(lft_m_) 0.015
## test_cndtn1 -0.041 0.005
## scl(l__):_1 0.003 -0.035 0.025# Take just the trials within condition
within_cond_pst_s2 <- s2_pst %>% filter(test_condition == "cognitive_cognitive" | test_condition == "overt_overt")
within_cond_pst_s2$test_condition <-
factor(within_cond_pst_s2$test_condition, levels=c("overt_overt", "cognitive_cognitive"))
contrasts(within_cond_pst_s2$test_condition) <- c(-.5, .5)
head(within_cond_pst_s2$test_condition)## [1] cognitive_cognitive cognitive_cognitive cognitive_cognitive
## [4] cognitive_cognitive cognitive_cognitive cognitive_cognitive
## attr(,"contrasts")
## [,1]
## overt_overt -0.5
## cognitive_cognitive 0.5
## Levels: overt_overt cognitive_cognitiveUse the comparable model
summary(m0_s2 <- glmer(resp_num ~ scale(left_min_right) + (0 + scale(left_min_right)|ID),
data=within_cond_pst_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: resp_num ~ scale(left_min_right) + (0 + scale(left_min_right) |
## ID)
## Data: within_cond_pst_s2
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7161.3 7181.7 -3577.7 7155.3 6616
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.1402 -0.7899 0.1255 0.7803 4.7614
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID scale(left_min_right) 1.507 1.228
## Number of obs: 6619, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.11126 0.02914 3.818 0.000134 ***
## scale(left_min_right) 1.35153 0.11345 11.912 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## scl(lft_m_) 0.014summary(m1_s2 <- glmer(resp_num ~ scale(left_min_right)*test_condition + (0 + scale(left_min_right) |ID),
data=within_cond_pst_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## resp_num ~ scale(left_min_right) * test_condition + (0 + scale(left_min_right) |
## ID)
## Data: within_cond_pst_s2
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7163.5 7197.5 -3576.8 7153.5 6614
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.0975 -0.7881 0.1258 0.7749 4.9593
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID scale(left_min_right) 1.509 1.228
## Number of obs: 6619, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.11101 0.02914 3.809 0.00014 ***
## scale(left_min_right) 1.35180 0.11350 11.910 < 2e-16 ***
## test_condition1 0.06205 0.05826 1.065 0.28689
## scale(left_min_right):test_condition1 -0.05321 0.06672 -0.797 0.42518
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sc(__) tst_c1
## scl(lft_m_) 0.014
## test_cndtn1 -0.008 0.004
## scl(l__):_1 0.006 -0.002 0.027p_s1 <-
sjPlot::plot_model(m1_s1, type="pred", terms = c("left_min_right [all]", "test_condition"))
p_fin_1 <- p_s1 + ap + tp + ga +
xlab("Left minus right reward history") + ylab("") + ggtitle("") +
theme(plot.title = element_text(hjust=0)) +
ylab("Chance of picking left") +
xlab("Left minus right reward history") +
ggtitle("Predicted probabilities") +
scale_color_manual(values=c("orange", "purple"), labels=c("overt-overt", "cognitive-cognitive")) + tol## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.p_s2 <-
sjPlot::plot_model(m1_s2, type="pred", terms = c("left_min_right [all]", "test_condition"))
p_fin_s2 <- p_s2 + ap + tp + ga +
xlab("Left minus right reward history") + ylab("") + ggtitle("") +
theme(plot.title = element_text(hjust=0)) +
ylab("Chance of picking left") +
xlab("Left minus right reward history") +
#ggtitle("Predicted probabilities") +
scale_color_manual(values=c("orange", "purple"), labels=c("overt-overt", "cognitive-cognitive")) +
theme(legend.text = element_text(size = 18),
legend.title = element_blank(),
legend.key.size = unit(1.3, 'lines'))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.p_fin_1 + p_fin_s2
Mixed test phase
unique(s1_pst$test_condition)## [1] "overt_cogn" "overt_overt" "cogn_cogn"mix_s1 <- s1_pst %>% filter(test_condition == "overt_cogn")# If cognitive is on the left and they chose left then they chose cognitive
mix_s1[mix_s1$left_training_cond=="cognTraining", "chose_cognitive"] <-
if_else(mix_s1[mix_s1$left_training_cond=="cognTraining", "response"]=="left", 1, 0)
# If cognitive is on the right and they chose right then they chose cognitive
mix_s1[mix_s1$right_training_cond=="cognTraining", "chose_cognitive"] <-
if_else(mix_s1[mix_s1$right_training_cond=="cognTraining", "response"]=="right", 1, 0)table(mix_s1$chose_cognitive)[2]/sum(table(mix_s1$chose_cognitive))## 1
## 0.4751506unique(s2_pst$test_condition)## [1] "overt_cognitive" "cognitive_cognitive" "overt_overt"mix_s2 <- s2_pst %>% filter(test_condition == "overt_cognitive")# If cognitive is on the left and they chose left then they chose cognitive
mix_s2[mix_s2$left_training_cond=="cognitive", "chose_cognitive"] <-
if_else(mix_s2[mix_s2$left_training_cond=="cognitive", "response"]=="left", 1, 0)
# If cognitive is on the right and they chose right then they chose cognitive
mix_s2[mix_s2$right_training_cond=="cognitive", "chose_cognitive"] <-
if_else(mix_s2[mix_s2$right_training_cond=="cognitive", "response"]=="right", 1, 0)table(mix_s2$chose_cognitive)[2]/sum(table(mix_s2$chose_cognitive))## 1
## 0.4231989mixed_summs_s1 <- mix_s1 %>% group_by(choice_type_plotting) %>%
summarize(m=mean(chose_cognitive), n())
mixed_summs_s1## # A tibble: 4 × 3
## choice_type_plotting m `n()`
## <chr> <dbl> <int>
## 1 "freq-P\nfreq-P" 0.490 997
## 2 "freq-R\nfreq-R" 0.428 997
## 3 "infreq-P\ninfreq-P" 0.553 992
## 4 "infreq-R\ninfreq-R" 0.429 998mixed_summs_s2 <- mix_s2 %>% group_by(choice_type_plotting) %>%
summarize(m=mean(chose_cognitive), n())
mixed_summs_s2## # A tibble: 4 × 3
## choice_type_plotting m `n()`
## <chr> <dbl> <int>
## 1 "freq-P\nfreq-P" 0.408 1104
## 2 "freq-R\nfreq-R" 0.393 1104
## 3 "infreq-P\ninfreq-P" 0.447 1102
## 4 "infreq-R\ninfreq-R" 0.445 1104summary(choice_effect_s1 <- glmer(chose_cognitive ~ 1 + ( 1 |ID),
data=mix_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: chose_cognitive ~ 1 + (1 | ID)
## Data: mix_s1
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 5269.0 5281.6 -2632.5 5265.0 3982
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4732 -0.8937 -0.4973 0.9600 2.3488
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.5492 0.7411
## Number of obs: 3984, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.10926 0.07443 -1.468 0.142summary(choice_effect_s2 <- glmer(chose_cognitive ~ 1 + ( 1|ID),
data=mix_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: chose_cognitive ~ 1 + (1 | ID)
## Data: mix_s2
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 5829.4 5842.2 -2912.7 5825.4 4412
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5760 -0.7917 -0.5754 1.0417 2.2498
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.3916 0.6258
## Number of obs: 4414, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.33769 0.06215 -5.433 5.53e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Questionnaires
s1_qs <- read.csv("../data/cleaned_files/s1_qs_deidentified.csv")
s2_qs <- read.csv("../data/cleaned_files/s2_qs_deidentified.csv")Exclude pts who answered both catch questions incorrectly
S1 all correct so no exclusions
table(s1_qs$Q6_20) # All correct##
## A little bit
## 123table(s1_qs$Q7_19) # All correct##
## Very false for me
## 123table(s2_qs$Q6_20) # All correct##
## A little bit
## 137table(s2_qs$Q7_19) # One pt incorrect but they answered the other catch item correctly, so retained (given recs in lit that excluding on a single catch item can bias results)##
## Very false for me Very true for me
## 136 1#s2_qs %>% filter(Q7_19 == "Very true for me")S1 questionnaires
Get and recode PTQ and RRS
s1_ptq <- s1_qs %>% select(contains("Q2_")) %>% setNames(paste0("PTQ_", 1:15))
s1_ptq$ID <- s1_qs$ID
s1_ptq_copy <- s1_ptq
s1_ptq_copy$ID <- s1_qs$ID
s1_ptq[s1_ptq=="Never"] <- 0
s1_ptq[s1_ptq=="Rarely"] <- 1
s1_ptq[s1_ptq=="Sometimes"] <- 2
s1_ptq[s1_ptq=="Often"] <- 3
s1_ptq[s1_ptq=="Almost always"] <- 4s1_rrs <-
s1_qs %>% select(contains("Q8_")) %>% setNames(paste0("RRS_", 1:10))
s1_rrs$ID <- s1_qs$ID
s1_rrs_copy <- s1_rrs
s1_rrs_copy$ID <- s1_qs$ID
s1_rrs[s1_rrs=="Almost never"] <- 1
s1_rrs[s1_rrs=="Sometimes"] <- 2
s1_rrs[s1_rrs=="Often"] <- 3
s1_rrs[s1_rrs=="Almost always"] <- 4Drop the 1 pt who skipped all items (leaving n=122 completers for study 1)
s1_ptq <- s1_ptq[-66, ]
s1_rrs <- s1_rrs[-66, ]Convert to numeric
s1_ptq_num <- foreach (i = 1:ncol(s1_ptq)) %do% {
data.frame(as.numeric(unlist(s1_ptq[i]))) %>% setNames(names(s1_ptq[i]))
}%>% bind_cols()Just 4 data points (< .1% of data) missing for remaining, so mean impute these
length(which(is.na(s1_ptq_num[, 1:15])))## [1] 4length(which(is.na(s1_ptq_num[, 1:15])))/(dim(s1_ptq_num[, 1:15])[1]*dim(s1_ptq_num[, 1:15])[2]) ## [1] 0.002185792s1_ptq_final <- foreach (i = 1:nrow(s1_ptq_num)) %do% {
if (any(is.na(s1_ptq_num[i, ]))) {
s1_ptq_num[i, is.na(s1_ptq_num[i, ])] <- mean(unlist(s1_ptq_num[i, ]), na.rm=TRUE)
}
s1_ptq_num[i, ]
}%>% bind_rows()hist(rowSums(s1_ptq_final[, 1:15]), breaks=25)
s1_ptq_final$ptq_sum <- rowSums(s1_ptq_final[, 1:15])Spot check IDs lined up properly
# s1_ptq_final %>% filter(ID == 22)
# s1_ptq_copy %>% filter(ID == 22)
# s1_qs %>% filter(ID == 22)
# s1_ptq_final %>% filter(ID == 104)
# s1_ptq_copy %>% filter(ID == 104)
# s1_qs %>% filter(ID == 104)Convert to numeric
s1_rrs_num <- foreach (i = 1:ncol(s1_rrs)) %do% {
data.frame(as.numeric(unlist(s1_rrs[i]))) %>% setNames(names(s1_rrs[i]))
}%>% bind_cols()length(which(is.na(s1_rrs_num[, 1:10])))## [1] 3length(which(is.na(s1_rrs_num[, 1:10])))/(dim(s1_rrs_num[, 1:10])[1]*dim(s1_rrs_num[, 1:10])[2]) ## [1] 0.002459016s1_rrs_final <- foreach (i = 1:nrow(s1_rrs_num)) %do% {
if (any(is.na(s1_rrs_num[i, ]))) {
s1_rrs_num[i, is.na(s1_rrs_num[i, ])] <- mean(unlist(s1_rrs_num[i, ]), na.rm=TRUE)
}
s1_rrs_num[i, ]
}%>% bind_rows()Spot check
# s1_rrs_final %>% filter(ID == 106)
# s1_rrs_copy %>% filter(ID == 106)
# s1_qs %>% filter(ID == 106)hist(rowSums(s1_rrs_final[, 1:10]), breaks=25)
s1_rrs_final$rrs_sum <- rowSums(s1_rrs_final[, 1:10])Sanity check RRS and PTQ are strongly correlated
ComparePars(s1_rrs_final$rrs_sum, s1_ptq_final$ptq_sum, use_identity_line = 0)
Reduced to just matching IDs
m1_s1_model_red <- m1_study1_eb %>% filter(ID %in% s1_ptq_final$ID)
m1_s1_model_red$c_min_o_phi <- m1_s1_model_red$cog_phi - m1_s1_model_red$overt_phiassert(all(m1_s1_model_red$ID==s1_ptq_final$ID))
assert(all(m1_s1_model_red$ID==s1_rrs_final$ID))PTQ
ComparePars(m1_s1_model_red$c_min_o_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)
ComparePars(m1_s1_model_red$cog_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)
ComparePars(m1_s1_model_red$overt_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)
Same basic pattern for RRS
ComparePars(m1_s1_model_red$c_min_o_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)
# Would expect higher decay to positively correlate w more brooding — instead small and ns neg
ComparePars(m1_s1_model_red$cog_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)
ComparePars(m1_s1_model_red$overt_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)
s2 questionnaires
Get and recode PTQ and RRS
s2_ptq <- s2_qs %>% select(contains("Q2_")) %>% setNames(paste0("PTQ_", 1:15))
s2_ptq$ID <- s2_qs$ID
s2_ptq_copy <- s2_ptq
s2_ptq_copy$ID <- s2_qs$ID
s2_ptq[s2_ptq=="Never"] <- 0
s2_ptq[s2_ptq=="Rarely"] <- 1
s2_ptq[s2_ptq=="Sometimes"] <- 2
s2_ptq[s2_ptq=="Often"] <- 3
s2_ptq[s2_ptq=="Almost always"] <- 4s2_rrs <-
s2_qs %>% select(contains("Q8_")) %>% setNames(paste0("RRS_", 1:10))
s2_rrs$ID <- s2_qs$ID
s2_rrs_copy <- s2_rrs
s2_rrs_copy$ID <- s2_qs$ID
s2_rrs[s2_rrs=="Almost never"] <- 1
s2_rrs[s2_rrs=="Sometimes"] <- 2
s2_rrs[s2_rrs=="Often"] <- 3
s2_rrs[s2_rrs=="Almost always"] <- 4No pts with skips
Convert to numeric
s2_ptq_num <- foreach (i = 1:ncol(s2_ptq)) %do% {
data.frame(as.numeric(unlist(s2_ptq[i]))) %>% setNames(names(s2_ptq[i]))
}%>% bind_cols()length(which(is.na(s2_ptq_num[, 1:15])))## [1] 3length(which(is.na(s2_ptq_num[, 1:15])))/(dim(s2_ptq_num[, 1:15])[1]*dim(s2_ptq_num[, 1:15])[2]) ## [1] 0.001459854s2_ptq_final <- foreach (i = 1:nrow(s2_ptq_num)) %do% {
if (any(is.na(s2_ptq_num[i, ]))) {
s2_ptq_num[i, is.na(s2_ptq_num[i, ])] <- mean(unlist(s2_ptq_num[i, ]), na.rm=TRUE)
}
s2_ptq_num[i, ]
}%>% bind_rows()hist(rowSums(s2_ptq_final[, 1:15]), breaks=25)
s2_ptq_final$ptq_sum <- rowSums(s2_ptq_final[, 1:15])Spot check IDs lined up properly
# s2_ptq_final %>% filter(ID == 95)
# s2_ptq_copy %>% filter(ID == 95)
# s2_qs %>% filter(ID == 95)Convert to numeric
s2_rrs_num <- foreach (i = 1:ncol(s2_rrs)) %do% {
data.frame(as.numeric(unlist(s2_rrs[i]))) %>% setNames(names(s2_rrs[i]))
}%>% bind_cols()length(which(is.na(s2_rrs_num[, 1:10])))## [1] 3length(which(is.na(s2_rrs_num[, 1:10])))/(dim(s2_rrs_num[, 1:10])[1]*dim(s2_rrs_num[, 1:10])[2]) ## [1] 0.002189781s2_rrs_final <- foreach (i = 1:nrow(s2_rrs_num)) %do% {
if (any(is.na(s2_rrs_num[i, ]))) {
s2_rrs_num[i, is.na(s2_rrs_num[i, ])] <- mean(unlist(s2_rrs_num[i, ]), na.rm=TRUE)
}
s2_rrs_num[i, ]
}%>% bind_rows()Spot check
# s2_rrs_final %>% filter(ID == 128)
# s2_rrs_copy %>% filter(ID == 128)
# s2_qs %>% filter(ID == 128)hist(rowSums(s2_rrs_final[, 1:10]), breaks=25)
s2_rrs_final$rrs_sum <- rowSums(s2_rrs_final[, 1:10])Sanity check RRS and PTQ are strongly correlated
ComparePars(s2_rrs_final$rrs_sum, s2_ptq_final$ptq_sum, use_identity_line = 0)
Reduced to just matching IDs
m1_s2_model_red <- m1_study2_eb %>% filter(ID %in% s2_ptq_final$ID)
m1_s2_model_red$c_min_o_phi <- m1_s2_model_red$cog_phi - m1_s2_model_red$overt_phiassert(all(m1_s2_model_red$ID==s2_ptq_final$ID))
assert(all(m1_s2_model_red$ID==s2_rrs_final$ID))PTQ
ComparePars(m1_s2_model_red$c_min_o_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)
ComparePars(m1_s2_model_red$cog_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)
ComparePars(m1_s2_model_red$overt_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)
Same basic pattern for RRS
ComparePars(m1_s2_model_red$c_min_o_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)
# Would expect higher decay to positively correlate w more brooding — instead small and ns neg
ComparePars(m1_s2_model_red$cog_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)
ComparePars(m1_s2_model_red$overt_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)